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Abstract 


Some  important  advances  took  place  during  the  last  several  years 
in  the  development  of  genuinely  multidimensional  upwind  schemes 
for  the  compressible  Euler  equations.  In  particular,  a  robust,  high- 
resolution  genuinely  multidimensional  scheme  which  can  be  used 
for  any  of  the  flow  regimes  computations  was  constructed  (see  [1- 
3]).  This  paper  summarizes  briefly  these  developments  and  outlines 
the  fundamental  advantages  of  this  approach. 
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1  INTRODUCTION 


The  efficiency  of  the  existing  steady-state  multigrid  solvers  routinely  xised  for 
the  flow  problems  in  engineering  practice  is  still  very  poor.  There  clearly  exists 
a  pressing  need  for  more  efficient  algorithms. 

In  order  to  obtain  a  truly  efficient  steady-state  solver,  some  fundamental  issues 
concerning  different  a^ects  of  an  algorithm  need  to  be  addressed.  The  recently 
proposed  genuinely  multidimensional  approach  towards  the  construction  of 
the  discrete  schemes  for  the  compressible  flow  resolves  some  of  these  issues.  A 
discussion  in  this  regard  is  the  main  .subject  of  this  paper. 


1.1  Genuinely  multidimensional  schemes 


The  que.st  for  a  genuinely  multidimensional  upwind  scheme  began  more  than 
a  decade  ago.  Initially  it  was  motivated  chiefly  by  the  esipectation  that,  once 
such  a  scheme  is  developed,  it  will  imitate  the  physics  of  the  fluid  flow  more 
accurately  than  the  standard  dimension-by-dimension  approach.  It  was,  how¬ 
ever,  suggested  later  in  [4], [5]  that  improving  the  efficiency  of  the  the  steady- 
state  .solvers  may  be  the  most  important  reason  for  developing  the  genuinely 
multidimensional  approach. 

One  of  the  Tnfl.in  difficulties  in  the  numerical  treatment  of  compressible  flow  is 
the  possible  pre.sence  of  .shocks  in  the  solution.  It  is  well  known  that  a  scheme 
that  is  both  .second  order  accurate  and  avoids  under-  and  overshoots  (which 
may  trigger  a  nonlinear  instability)  near  discontinuities  has  to  be  nonlinear. 
Such  a  scheme  has  to  incorporate  the  so-called  high-resolution  mechanism, 
i.e.  a  smoothness  monitor,  that  is  usually  implemented  in  the  form  of  a  flux- 
limiter.  Initially,  such  schemes  were  developed  for  the  one-dimen.sional  case. 
Then,  extending  this  approach  to  multidimensions  was  done  on  a  dimension- 
by-dimension  basis.  The  well  known  fact,  however,  is  that  the  Gauss-Seidel  re¬ 
laxation  is  unstable  when  applied  in  conjunction  with  such  schemes  [6] .  There¬ 
fore,  the  standard  multigrid  solvers  have  to  re.sort  to  the  defect-correction 
technique  or  multistage  Runge-Kutta  relaxation  and  the  efficiency  of  such 
.solvers  may  be  poor.  A  closer  look  reveals  that  the  .standard  high-resolution 
discretizations  suffer  from  the  following  deficiency:  the  high-frequency  error 
components  may  be  (nearly)  invisible  to  the  residuals  of  the  discrete  equa¬ 
tions,  i.e.  the  discrete  scheme  is  (nearly)  unstable.  In  turn  this  means  that 
it  may  be  inherently  impossible  to  construct  a  good  .smoother  (an  important 
ingredient  of  a  multigrid  solver)  using  these  discrete  schemes. 

A  genuinely  multidimen.sional  advection  scheme  was  constnicted  in  [4,5].  The 
scheme  was  named  “genuinely  multidimensional”  since  it  imitates  well  the 
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anisotropy  of  the  advection  phenomena  in  two  dimensions:  artificial  dissipation 
is  added  only  along  the  streamline,  while  the  high-resolution  mechanism  affects 
significantly  the  aoss-fiow  direction  only.  The  key  feature  of  this  scheme  is 
the  two-dimensional  limiter,  i.e.  the  argument  of  a  limiter-function  is  the  ratio 
of  finite  differences  in  two  different  coordinate  directions.  The  scheme  was 
formulated  in  the  control-volume  context  for  Cartesian  grids  and  relied  on 
the  compact  9-point-box  .stencil.  The  fundamental  advantage  of  this  approach 
is  that  the  two-dimensional  high-resolution  mechanism  does  not  damage  the 
stability  properties  of  the  discretization. 

The  so-called  “residual  distribution”  (or  “fluctuation-splitting”)  schemes  for 
scalar  advection  equation  on  unstructured  triangular  grids  were  presented  in 
[7].  It  was  found  later  that  the.se  schemes  have  some  links  to  the  aforemen¬ 
tioned  genuinely  two-dimensional  control-volume  approach  for  the  advection 
equation.  Exploration  of  these  links  led  to  the  unification  of  the  two  approaches 
and  re.sulted  in  a  scheme  that  incorporated  two-dimen.sional  limiters  and  was 
formulated  for  unstnictured  triangular  grids.  This  scheme  (like  that  presented 
in  [4], [5])  can  be  given  a  purely  algebraic  interpretation.  However,  the  task  of 
extending  the.se  ideas  to  systems  of  equations  appeared  to  be  a  complicated 
one. 

Consider  a  hjq)erbolic  system  of  partial  differential  equations  in  two  dimen¬ 
sions 


Ut  +  AUj;  +  BUy  =  0, 


(1) 


where  u  is  the  vector  of  size  N  and  A,  B  are  N  x  N  matrices.  The  matrices 
A  and  B  in  general  do  not  commute.  This  means  that  they  cannot  be  di¬ 
agonalized  simultaneously,  i.e.  the  sy.stem  cannot  be  written  as  N  advection 
equations. 

A  prolonged  effort  was  to  represent  locally  the  physics  of  compressible  flow  by 
finite  number  of  simple  waves  using  the  tocal  gradients  of  the  solution  (in  the 
spirit  of  [8])  with  intention  to  apply  a  genuinely  two-dimensional  advection 
scheme  to  each  one  of  the  simple  waves.  However,  the  schemes  constructed  in 
this  way  for  the  Euler  equations  suffered  from  a  severe  lack  of  robustness. 

The  breakthrough  approach  that  resulted  in  a  robust  genuinely  multidimen¬ 
sional  scheme,  suitable  for  the  computations  of  the  entire  range  of  the  flow 
regimes  was  pre.sented  in  [1] .  Then  it  was  described  in  more  detail  including  the 
implications  for  multigrid  and  extension  to  3D  in  [2]  and  [3].  The  key  idea  was 
not  to  try  to  apply  the  multidimensional  advection  schemes  to  systems,  but 
rather  the  same  strategy  that  was  u.sed  to  construct  the  scalar  scheme.  The  al¬ 
gebraic  interpretation  of  the  advection  scheme  played  an  important  role  at  this 
point.  It  was  crucial  to  recognize  that  a  certain  linear  first  order  scheme  based 
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on  standard  upwind  methodology  can  be  used  as  a  basic  building  block  for  the 
hyperbolic  systems  as  well  as  for  the  scalar  advection.  The  multidimensional 
high-resolution  corrections  are  then  applied  in  a  formal  way  similarly  to  the 
scalar  case.  The  resulting  scheme  for  the  Euler  equations  was  demonstrated 
to  produce  a  very  good  quality  solution  for  subsonic,  transonic  and  s\ipersonic 
regimes.  The  approach  was  called  “genuinely-mtiltidimensional”  since  it  can 
be  argued  that  it  leads  to  a  discrete  scheme  whose  artificial  dissipation  is  a 
rotationaHy-invariant  differential  operator  {in  other  words  -  the  artificial  dis¬ 
sipation  operator  is  independent  to  a  certain  extent  of  the  grid  direction).  It 
was  not  clear  if  this  particular  property  is  of  any  direct  practical  importance. 

The  constructed  high-resolution  scheme  for  the  Euler  equations  relies  on  a 
compact  stencil.  The  result  of  this  property  is  a  smaller  error  in  smooth  regions 
and  better  resolution  of  discontinuities,  comparing  to  the  standard  dimension- 
by-dimension  approach.  However,  in  our  view,  these  are  only  marginal  im¬ 
provements.  The  fundamental  advantage  of  this  approach  is  that  it  leads  to 
a  scheme  that  combines  high-resolution  and  good  stability  properties.  It  was 
demonstrated  in  [2,3]  that  the  Collective  Gauss-Seidel  relaxation  is  stable 
when  applied  directly  to  the  resulting  high-resolution  discretization  of  the  hy¬ 
perbolic  systems.  This  results  in  a  very  simple,  efiicient  and  robust  multigrid 
solver  for  the  compressible  Euler  equations,  suitable  for  the  entire  range  of 
flow  regimes. 

Some  researchers  who  were  previoiisly  pursuing  other  directions  adopted  the 
genuinely-multidimensional  approach  proposed  in  [1—3]  and  attributed  to  it 
a  term  “Positive  Matrix  Distribution  Schemes”.  A  modification  of  the  un¬ 
derlying  first  order  scheme  (the  system  N-scheme)  aimed  at  improving  the 
discontinuity  resolution  was  proposed  by  van  der  Weide  and  Deconinck  in  [9]. 

It  should  be  mentioned  that  important  steps  towards  the  construction  of  a 
genuinely  mxiltidimensional  schemes  for  the  Euler  equations  were  made  by 
Colella  [10,11],  LeVeque  [12]  and  Radvogin  [13].  However,  the  nonlinear  high- 
resolution  corrections  in  these  schemes  rely  on  one-dimensional  limiters,  which 
introduces  some  of  the  dimension-by-dimension  flavor. 

Another  very  interesting  approach  was  proposed  in  [14].  A  discretization  for 
the  triangular  unstructured  grids  for  the  Euler  equations  was  developed.  The 
problem  was  that  the  scheme  was  linear,  i.e.  it  did  not  incorporate  any  non¬ 
linear  high-resolution  mechanism. 


1.2  Multigrid  for  advection  dominated  problems 


One  of  the  major  reasons  for  the  poor  efficiency  of  the  standard  flow  solvers 
(see  [15])  is  the  fact  that  for  advection  dominated  problems  the  coarse  grid 
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provides  only  a  fraction  of  the  needed  correction  for  certain  error  components. 
It  is  well  known  that  the  steady  Euler  equations  contain  two  different  factors: 
the  advection  and  the  Full-Potential  type  operators.  The  latter  is  either  of 
elliptic  or  hyperbolic  type  depending  on  the  flow  r^me  (subsonic  or  super¬ 
sonic).  The  difficulty  mentioned  above  can  be  avoided  ([15])  by  constructing 
a  solver  that  distinguishes  between  different  factors  of  the  system  and  treats 
each  one  appropriately.  In  the  subsonic  case,  for  instance,  the  advection  factor 
can  be  treated  by  marching  and  the  elliptic  factor  by  multigrid.  The  efficiency 
of  such  an  algorithm  will  be  essentially  the  same  eis  that  of  the  multigrid 
solver  for  the  elliptic  part  only.  Such  algorithms  are  referred  to  as  “essentially 
optimal” .  An  approach  to  achieve  a  separation  between  the  co-factors  -  the 
so-called  Distributive  Gauss-Seidel  relaxation  -  was  proposed  in  [15].  It  was 
demonstrated  in  [16]  that  using  this  approach  one  can  obtain  the  essentially 
optimal  multigrid  efficiency  for  a  staggered-grid  discretization  of  the  incom¬ 
pressible  Navier-Stokes  equations.  It  is  interesting  to  mention  that  a  similar 
observation  was  made  earlier  in  [17]. 

A  related  approach  was  proposed  by  Ta’asan  [18]  for  the  incompressible  and 
compressible  subsonic  Euler  equations.  The  staggered-grid  discretization  is 
based  on  the  canonical  variables  formulation  (see  [19]),  that  expresses  the 
partitioning  of  the  steady  Euler  equation  into  elliptic  and  advection  factors. 
The  essentially  optimal  multigrid  efficiency  was  demonstrated  using  this  ap¬ 
proach  for  subsonic  flow  and  body-fitted  grids.  A  possible  limitation  of  this 
approach  may  be  that  it  is  not  directly  generalizable  for  the  viscous  flow. 

Another  way  towards  achieving  the  optimal  multigrid  efficiency  is  based  on  the 
pressure  equation  formulation  of  the  Euler  equations.  We  describe  it  briefly  in 
this  paper  as  well.  This  approach  is  based  on  a  very  old  idea  and  is  a  general¬ 
ization  [20].  Its  main  virtue  is  simplicity.  It  can  also  be  classified  as  Weighted 
Gauss-Seidel  relaxation  [15].  An  extensive  set  of  numerical  compiitations  us¬ 
ing  this  scheme  together  with  more  details  regarding  the  implementaion  is 
reported  in  [21]).  The  limitation  of  this  approach,  however,  is  that  it  is  not 
clear  so  far  if  it  is  generalizable  to  viscous  compressible  flow. 

Following  the  work  of  Ta’asan,  some  researchers  attempted  to  apply  the  idea  of 
partitioning  the  Euler  equations  towards  the  construction  of  discrete  schemes 
(see  [22]  and  [23]).  It  is  well-known  that  the  two-dimensional  Euler  equations 
in  supersonic  case  can  be  written  as  four  locally  decoupled  advection  equations 
(see  [24]).  This  property  was  used  as  a  basis  for  applying  the  advection  schemes 
to  discretize  the  system  in  this  case.  In  siibsonic  case,  however,  the  distinc¬ 
tion  was  made  between  the  advection  and  the  elliptic  ( “acoustic  subsystem” ) 
partitions.  The  treatment  of  transonic  flow,  however,  was  problematic  since  it 
required  matching  of  two  different  discretizations  accross  the  sonic  lines.  An¬ 
other  drawback  of  these  approaches  was  that  they  are  cannot  be  generalized 
to  three  dimensions  (see  [9]).  To  conclude,  the  discrete  schemes  constructed  in 


4 


this  way  suffered  from  a  lack  of  robustness  and  generality.  No  optimal  multi- 
grid  eflBciency  was  demonstrated  either.  Finally,  some  of  the  researchers  who 
previously  followed  this  direction  adopted  the  genuinely-multidimensional  ap¬ 
proach  proposed  in  [1-3]  (.see  [9]). 


1.3  What  this  paper  is  about 


In  this  paper  we  first  present  a  brief  review  of  the  construction  of  genuinely 
multidimensional  schemes  for  the  scalar  advection  and  the  compressible  Euler 
equations.  We  siimmarize  the  basic  properties  of  the  discretizations,  emphasiz¬ 
ing  those  that  are  unique  to  this  approach  and  are  of  fundamental  importance 
for  practical  purposes. 


The  separation  of  the  co-factors  can  be  in  general  achieved  in  two  ways.  One 
way  is  to  cast  equations  into  such  a  form  that  the  different  co-factors  can  be 
discretized  separately.  The  canonical  variables  approach  by  Ta’asan  [19]  can 
be  ela-ssified  as  such.  The  pressure  equation  based  scheme,  presented  briefiy  in 
this  paper,  belongs  to  this  category  as  well.  The  key  advantage  of  this  type  of 
methods  is  in  their  simplicity  (this  is  e.specially  tnie  for  the  pressure  equation 
based  scheme).  The  disadvantage,  though,  may  be  in  the  unsuflacient  gener¬ 
ality.  Another  more  general  way  to  achieve  the  optimal  multigrid  efficiency 
(see  [15])  is  to  discretize  the  equations  in  .some  primitive  form  and  to  apply 
a  relaxation  of  the  Distrib\itive  Gau.s.s- Seidel  type.  Such  relaxation  .should  be 
designed  in  such  a  way  that  it  distinguishes  between  the  different  co-factors 
of  the  system  and  treats  each  one  of  them  appropriately.  It  appears,  however, 
that  in  order  to  achieve  this  not  any  discrete  schemes  are  .suitable,  but  only 
those  satisfying  a  certain  condition.  As  it  was  mentioned  before,  it  is  not  clear, 
what  are  the  direct  practical  implications  of  the  genuine  multidimensionality 
(or  the  rotational  invariance  of  the  artificial  dissipation)  property  of  the  ap¬ 
proach.  However,  we  present  in  this  paper  a  heuristc  argument  suggesting  that 
the  genuine  multidimensionality  is  closely  related  to  the  factorizability  prop¬ 
erty  of  the  discretization.  The  latter  is  of  fundamental  practical  importance. 
It  is  necessary  in  order  to  construct  a  Distributive  Gau.ss-Seidel  relaxation 
that  will  allow  to  decouple  the  advection  and  Full-Potential  co-factors  of  the 
Euler  system  and  thus  to  obtain  an  optimally  efficient  multigrid  solver.  This 
approach  is  very  general  since  it  does  not  require  to  cast  the  equations  into 
any  special  form. 
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Fig.  1.  Triangle. 

2  MULTIDIMENSIONAL  UPWINDING 

In  this  section  we  briefly  review  the  construction  of  the  high-resolution  gen¬ 
uinely  multidimensional  upwind  schemes  for  the  scalar  advection  equation  and 
for  the  Euler  system.  Consider  a  general  triangular  grid  covering  the  domain. 
Assume  the  discrete  unknowns  are  deflned  at  the  grid  nodes  thus  defining  a 
linear  function  and  allowing  us  to  evaluate  the  gradients  of  the  current  so¬ 
lution  approximation  on  each  triangle.  The  approach  is  to  construct  discrete 
equations  which  are  to  be  solved  at  each  node  from  the  portions  of  the  residu¬ 
als  of  the  equations  computed  on  the  triangles  having  this  node  as  a  common 
vertex.  In  other  words,  portions  of  the  residual  of  the  equation  evaluated  on  a 
particular  triangle  contribute  to  the  construction  of  the  discrete  equations  at 
the  nodes  of  this  triangle.  The  problem  is  to  find  the  exact  rules  for  this  con- 
stnrction,  so  that  the  resulting  discrete  equations  will  have  certain  desirable 
properties. 


2.1  Advection  scheme 


We  consider  triangular  element  T,  and  choose  two  out  of  the  three  faces.  Then 
we  write  the  advection  equation  we  wish  to  solve  in  the  local  coordinate  system 
aligned  with  these  two  faces 

Ut  -h  aux  +  buy  =  0  (2) 


Without  loss  of  generality  we  can  consider  a  linear  constant  coefficient  equa¬ 
tion,  since  in  general  non-linear  case  we  can  linearize  the  equation  on  each 
triangle  ([25]). 

We  can  write  the  discrete  equation  at  the  grid  node  i  in  the  following  form 

^  (^) 

j^i 
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Definition  1  The  discrete  scheme  (3)  is  said  to  be  of  the  positive  type  if  Cij  > 
0  for  all  j . 

Numerical  soltuions  obtained  using  a  positive  scheme  satisfy  a  discrete  max- 
imimum  principle.  This  property  is  useful  to  ensure  that  the  discrete  solution 
will  be  non-oscillatory  near  discontinuities. 

We  shall  oTitline  here  the  construction  of  a  positive  advection  scheme.  Residual 
of  the  eqtiation  (2)  can  be  represented  as  a  sum  of  two  portions 

r  =  r®  +  (4) 

where 

r®  =  -Srau^]  r^  =  Srlni^.  (5) 


Residual  of  the  equation  (2)  on  the  triangle  T  contributes  to  the  constniction 
of  the  discrete  equations  to  be  solved  at  each  of  the  three  nodes  of  the  triangle 
according  to  the  following  residual  distribution  formulae 

node  1  <—  r®(l  -  sign(a)) 

node  2  <—  r®(l  +  sign(a))  +  r*'(l  —  sign(6))  (6) 

node  3  <—  r^(l  +  sign(6)) 


It  easy  to  see  that  this  constniction  results  in  a  positive  scheme  since  for  any 
real  number  2  we  have  the  following  inequality  ±z{l  ±  sign(2r))  >  0.  The 
accuracy  of  such  a  scheme,  though,  is  only  first  order. 

Definition  2  A  discrete  scheme  is  called  linearity  preserving  if  whenever  the 
residual  r  on  the  triangle  T  vanishes,  the  contributions  due  to  this  residual 
lead  to  a  zero  update  of  the  solution  at  each  of  the  three  nodes  of  the  triangle. 

A  linearity  preserving  scheme  is  second  order  accurate. 

Define  the  following  quantities 

=  r®  +  r'^^(q)]  =  r^  +  r^^(q)/q  (7) 

where  q  =  — r®/r".  In  this  paper  we  assume  that  ^  is  the  minmod  limiter.  Sub¬ 
stituting  r®* ,  r^*  instead  of  r®,  r*'  into  (6)  we  obtain  a  high-resolution  scheme. 
The  important  feature  here  is  the  two-dimensionality  of  the  limiter,  i.e.  the  fact 
that  the  argument  of  the  limiter-function  is  a  ratio  of  numerical  derivatives  in 
two  different  coordinate  direction  ([5], [4]). 
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(8) 


Using  the  following  limiter  identity 

=  -r^^{q)/q, 


it  is  easy  to  verify  that  the  constructed  nonlinear  scheme  Ls  indeed  both  positive 
and  linearity  preserving. 


2. 2  Extension  to  the  Euler  system 


Consider  a  hyperbolic  system  of  partial  differential  equations 

Ut  +  Aux  +  Buy  =  0.  (9) 


The  discrete  equation  approximating  (9)  at  node  i  can  be  written  as  follows 

Ci,iUi  -  ^  CijUj  =  0,  (10) 

jjii 

Property  of  positivity  can  be  formally  extended  to  the  system  case. 

Definition  3  The  discrete  scheme  (10)  is  said  to  be  of  the  positive  type  if  the 
matrices  Ci^j  (for  all  j )  have  non-negative  eigenvalues. 

It  is  not  clear,  however,  how  to  generalize  the  maximum  principle  for  sys¬ 
tems.  No  conclusions  can  be  derived  from  the  positivity  property  unless  the 
additional  assumption  about  the  symmetry  of  the  matrices  Gij  is  made.  The 
energy  stability  property  of  the  scheme  can  be  demonstrated  in  this  case. 
However,  for  the  Euler  system,  this  would  reqihre  the  use  of  the  symmetriz¬ 
ing  variables  formulation,  which  is  non-conservative.  This  makes  the  energy 
stability  property  too  restrictive  to  be  of  substantial  practical  importance  for 
the  Euler  equations. 

It  is  interesting  to  note  though  that  the  standard  high-resolution  schemes,  if 
carefully  implemented,  are  of  the  positive  type.  Therefore,  we  aimed  at  con- 
stmcting  a  genuinely  multidimensional  high-resolution  scheme  ([1-3]),  which 
is  of  the  positive  type  as  well. 

Assume  that  the  hyperbolic  system  of  equations  (9)  is  written  in  the  non- 
orthogonal  coordinate  frame  aligned  with  the  two  of  the  faces  of  triangle  T 
(Fig.l).  Residuals  of  the  system  on  triangle  T  can  be  represented  as  a  sum  of 
two  portions 

R  =  R”  +  R^  (11) 
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where 


R"  =  -5't^u^; 


(12) 


Consider  the  following  residual  distribution  formula 
node  1  R®(/  —  sign(A)) 

node  2  <-  R®(/  +  sign(^))  +  R^(/  -  sign(S))  (13) 

node  3  ^  R^(/  +  sign(R)) 


Assuming  that  matrix  M  has  a  complete  set  of  real  eigenvalues  (definition  of 
the  hyperbolicity  of  a  system)  it  is  easy  to  see  that  matrix  ±M{I  ±  sign(M)) 
is  non-negative  definite.  This  means  that  the  scheme  defined  by  (13)  is  of  the 
positive  type  by  construction  (as  an  upwind  scheme  is  expected  to  be). 

In  order  to  obtain  a  positive  high-resolution  genuinely  multidimensional  scheme 
for  a  hyperbolic  system  we  may  have  first  to  rewrite  system  (9)  (as  it  was  done 
in  [1-3])  in  a  different  set  of  variables.  The  auxiliary  variables  v  =  (s,  u,  v,  pY 
(see  [1-3])  are  a  good  choice  for  the  Euler  system  (here  s  is  the  entropy,  u,v 
-  the  velocity  components  and  p  is  the  pressure).  System  (9)  rewritten  in 
varibales  v  takes  the  following  form 

vt  -f-  Av^  -f  BYy  =  0  (14) 


where 


u  =  Tv 


(15) 


As  before,  we  can  comptite  residual  r  of  system  (14)  on  triangle  T  and  represent 
it  as  a  sum  of  two  portions: 


r  =  r  -H  r" 


(16) 


where 


r="  =  -StAv^] 

ry  = 


(17) 
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Considering  the  following  residual  distribution  formula 


node  1  -f—  Tr®(/  —  sign(^)) 

node  2  <-  T[r®(/  +  sign(^))  +  r^(/  -  sign(;B))]  (18) 

node  3  e-  Tr*'(/  +  sign(jB)), 


we  airive  at  a  positive  first  order  accurate  scheme  that  is  ideMical  to  (13). 
Introduce  the  following  quantities 

rf  =  rf  +  rf^{7,);  rf  =  +  rf  ^(g,)/g,  (19) 

with  qi  =  -rf/rf,  and  i  =  l,...,iV  Denote  by  r®*  and  r^*  vectors  whose 
components  are  rf  and  rf  {i  =  I,..., N)  respectively.  The  high-resolution 
genuinely  two-dimensional  scheme  can  be  obtained  by  substituting  r®  and 
instead  of  the  r®  and  r®'  into  (18). 

Using  the  limiter  identity 

rm<li)  =  -rmqi)/qi,  for  i  =  1,  N  (20) 


it  is  possible  to  show  that  the  genuinely  two-dimensional  high-resolution  scheme 
is  both  positive  and  linearity  preserving.  We  emphasize  here  again,  that  in 
order  to  achieve  this  property  for  the  Euler  equations,  it  was  necessary  to 
use  another  (non-conservative)  form  of  the  equations  when  introducing  the 
high-resolution  mechanism  ([1-3]).  The  auxiliary  variables  formulation  was 
suitable  for  this  purpose.  The  only  justification  for  the  desirability  of  the  pos¬ 
itivity  property  for  systems  of  equations  is  that  the  standard  high-resolution 
schemes,  if  properly  implemented,  have  this  property. 

Remark  4  Recall  that  when  constructing  a  discrete  scheme  we  have  chosen 
arbitrarily  two  out  of  three  faces  of  the  trinagle  T  (see  Fig.  1 )  for  the  purpose  of 
numerical  derivatives  evaluation.  Note,  that  any  choice  will  result  in  a  scheme 
with  good  properties.  In  the  case  of  the  scalar  advection,  though,  if  we  use 
two  faces  which  are  both  of  the  same  kind  with  respect  to  the  flow  direction 
(inflow  or  outflow),  the  resulting  scheme  will  rely  on  the  narrowest  possible 
stencil.  Therefore,  it  will  have  a  smaller  cross-stream  error  coefficient  and  wUl 
provide  a  somewhat  better  disontinuity  resolution.  It  is  possible  in  the  case  of 
the  Euler  system  to  optimize  the  resolution  of  a  specific  sharp  layer  (shock, 
contact  discontinuity  or  shear  layer)  by  the  appropriate  choice  of  the  faces 
(two  “closest”  to  the  layer  direction.  In  general,  though,  it  seems  reasonable 
to  choose  two  faces  which  are  either  inflow  or  outflow.  This  will  provide  better 
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resolution  of  the  contact  discontinuities  and  shear  layer,  while  shocks  are  well 
resolved  anyway. 

In  addition  to  the  smaller  error  in  smooth  regions  and  sharper  resolution  of 
discontinuities  (compared  to  the  standard  dimension-by-dimension  methods), 
the  constructed  scheme  also  offers  a  possibility  to  optimize  the  stencil  in  order 
to  resolve  a  particular  discontinuity  layer  (by  choosing  two  faces  of  a  triangular 
element  that  are  the  closest  to  the  direction  of  this  layer).  The  important 
advantage  of  this  approach,  however,  is  that  the  genuinely  multidimensional 
approach  presented  here  leads  to  a  scheme  which  has  both  high-resolution  and 
good  stability  properties. 


3  SEPARATION  OF  CO-FACTORS 


Incompressible  steady  Euler  equations  in  two  dimensions  can  be  written  in 
the  following  matrix  form  (assuming  that  the  fluid  density  p  =  1) 


Lu  =  0, 


(21) 


where 


u  =  {u,v,pf, 


L  = 


^  Q  0 

0  Q  dy  , 


ydxdy  0  J 


(22) 


(23) 


and  Q  =  U  V  is  the  advection  operator.  In  order  to  find  out  what  is  the  type 
the  system  of  equations  we  can  look  at  the  determinant  of  matrix  L 

det{L)  =  -Q^A  (24) 


The  co-factors  of  the  determinant  are  of  the  advection  and  another  of  the 
elliptic  types. 

Compressible  steady  Euler  equations  in  two  dimensions  can  be  written  in  the 
following  matrix  form 


Lu  =  0, 


(25) 
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where 


u  =  {s,u,v,pf, 

0  0  0 
0  pQ  0  dx 
Q  ^  pQ  dy 

^  0  pdy  J 


(26) 


(27) 


The  determinant  of  matrix  L  in  this  case 


dei(L)=/?2(Q)2[iQ2_A] 


(28) 


The  determinant  has  two  distinct  co-factors;  one  of  advection  and  another  of 
the  Full-Potential  type. 

In  the  rest  of  this  paper  we  shall  consider  the  subsonic  regime  only.  The  Full- 
Potential  factor  in  this  case  is  of  the  elliptic  type. 

It  was  suggested  in  [15]  that  the  different  co-factors  should  also  be  treated 
differently  (each  one  in  the  appropriate  way).  One  way  to  do  this  is  to  cast  the 

equations  into  such  a  form  that  the  co-factors  can  be  discretized  separately. 
Canonical  variables  formulation  [18]  as  well  as  the  pressure  equation  based 
schemes,  that  are  described  in  the  next  section,  belong  to  this  category.  A  more 
general  approach  (as  proposed  by  Brandt  in  [15])  is  to  discretize  the  equations 
in  some  primitive  form,  but  to  design  relaxation  of  Distributive  Gauss-Seidel 
type  (DCS)  such  that  it  will  separate  the  treatment  of  co-factors.  However,  the 
discrete  scheme  suitable  for  this  purpose  shotild  be  factorizahle,  i.e.  satisfy  the 
discrete  analog  of  the  property  (28).  We  show  in  this  paper  that  the  genuinely 
multidimensional  upwind  approach  leads  to  a  factorizahle  scheme. 


4  PRESSURE  EQUATION  APPROACH 


In  this  section  we  briefly  describe  the  pressure  equation  based  approach  flrst 
for  incompressible  and  then  for  compressible  cases. 
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4-1  Incompressible  case 


Considering  a  triangular  (unstnictured)  grid  and  assuming  that  the  unknowns 
are  located  at  grid  vertices,  we  can  define  piecewize-linear  functions  approxi¬ 
mating  each  tmknown.  Residual  of  the  Euler  equations  can  then  be  evaluated 
on  each  triangular  element  of  the  grid: 

r  =  L'‘u^  (29) 


where 


0 

0 

0  j 


and 

Q^  =  u^d^+v^d^ 


(30) 


(31) 


is  the  discrete  advection  operator. 

We  would  like  then  to  construct  the  discrete  equations  to  be  solved  at  a  certain 
grid  node  by  assembling  in  a  certain  way  the  residuals  of  the  Euler  equations 
comptited  on  the  triangular  element  having  this  node  as  a  common  vertex. 
These  equations  can  be  written  in  the  following  form 

Ph  =  P^L’^u'^  =  0  (32) 


In  order  to  obtain  the  equations  with  desirable  properties,  i.e.  to  achieve  the 
decoupling  of  the  co-factors,  we  can  define  the  following  assembly  matrix  P^ 


ph 


0/0 


(33) 


where 


(34) 
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is  the  operator  adjoint  to  Q^.  Then 


+  s.p.t. 


(35) 


where  is  the  discrete  Laplacian  and  “s.p.t.”  stands  for  snbprincipal  terms 
which  appear  due  to  the  non-constancy  of  the  advection  coefficients  and  are 
not  important  for  the  purpose  of  constructing  a  relaxation  scheme  [15].  The 
maxtrix  of  the  finite  difference  operators  in  (35)  is  upper-triangular.  The  pres¬ 
sure  is  subject  to  a  Poisson  equation  wich  is  decoiipled  (upto  subprincipal 
terms)  from  the  rest  of  the  system.  The  standard  Gaiiss-Seidel  relaxation 
can,  therefore,  be  applied  to  it  as  a  smoother.  The  momentum  equations  can 
be  looked  at  as  advection  equations  with  (known)  forcing  terms  (pressure 
derivatives).  The  ideas  of  multidimensional  upwinding  can  be  applied  to  dis¬ 
cretize  them.  The  strategy  applied  to  relax  the  system  can  be  then  to  relax 
the  pressure  first  and  then  to  update  the  velocities  components  by  relaxing 
the  momentum  equations  in  the  downstream  direction.  An  extensive  set  of 
computational  results  using  this  approach  is  presented  in  [21]. 


4..2  Compressible  case 


Similarly  to  the  incompressible  case,  residuals  of  the  equations  can  be  evalu¬ 
ated  on  each  triangular  element: 

r  =  L^u^,  (36) 


where  =  {s^,u^,v^,p^Y  and 

( a  0  0  ^ 

0  0 
0  0 
^  0  p^d^ 


(37) 


and  is  the  (discrete)  speed  of  sound.  Again,  we  constnict  the  discrete  equa¬ 
tions  to  be  solved  at  each  node  assembling  in  a  certain  way  the  residuals  of 
the  equations  on  the  elements  surrounding  this  node 

Ph  =  =  0  (38) 
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Choosing  the  assembly  matrix  to  be 


ph  _ 


^100  0  ^ 
0/0  0 
0  0/  0 
yOd^d^iQ*)^) 


(39) 


we  obtain 


Ph  = 


0  0 

0  /g'^  0 

0  0  p^Q’^ 

0  0  0 


0  ^ 

(s^\ 

yh 

+  s.p.t. 


(40) 


The  principal  part  of  last  equation  is  the  discrete  Prandtl-Glauert  (or  Full 
Potential)  operator  acting  on  the  pressxire.  This  operator  is  of  the  elliptic  type 
in  the  stibsonic  regime.  The  solution  process  of  the  resulting  discrete  equa¬ 
tions  is  very  similar  to  that  for  the  incompressible  case.  Some  numerical  tests 
illustrating  the  efficiency  of  the  multigrid  solver  based  on  this  discretization 
are  presented  in  the  §6.2. 


5  UPWINDING  AND  CO-FACTORS  SEPARATION 


Now  we  return  to  the  genuinely  multidimensional  upwind  approach.  We  would 
like  to  construct  a  linear  first  order  upwind  “positive”  scheme  such  that  it  is 
factorizable  and  is  also  upgradable  to  second  order  using  the  genuinely  two- 
dimensional  high-resolution  mechanism. 

First,  we  shall  take  a  closer  look  at  the  one-dimensional  case. 


5.1  One- dimensional  case 


Consider  a  first  order  upwind  scheme  for  the  one-dimensional  Euler  equations. 
Without  loss  of  generality  we  consider  the  primitive  variable  formulation 

=  0  (41) 
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where 


(42) 


and 


V  0 


+  g^'^)  df  -  , 

p{df-¥A^)  + 


(43) 


where  h  is  a  meshsize,  is  a  central  approximation  of  the  second  derivative, 
is  a  central  approximation  of  the  first  derivative  and  =  udl^  is  the 
advection  operator. 


5.1.1  Factorization 
The  determinant  of  L^: 

det{L^)  -  -  M^)d^J  (44) 


The  first  factor  is  the  upwind  scheme  approximating  an  advection  operator 
corresponding  to  the  entropy  equations.  The  Full-Potential  factor  is  approxi¬ 
mated  by  a  “short”  central  difference.  The  issue  of  factorization  appears  to  be 
trivial  in  this  case,  since  the  momentum  and  the  pressure  equations  correspond 
solely  to  the  elliptic  factor. 


5.1.2  Distributive  Gauss-Seidel  relaxation 
Introducing  new  variables 

(s'^,  u^, =  =  =  M'‘(s'‘,  4>^f  (45) 

the  Euler  system  will  take  the  following  form 

=  0  (46) 
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i+1 

Fig.  2.  One-dimensional  grid 

Assuming  that  that 


= 


/l  0  0  \ 

Vo  -p(af  -  Pilcdl  -  Q^^)J 


(47) 


we  obtain 


= 


/ -h\u\di^  +  Q^’^  0  0 

0  p{l-M^)d^,  0 

0  0 


(48) 


The  philosophy  of  the  Distributive  Gauss-Seidel  relaxation  is  as  follows:  We 
would  like  to  store  at  the  gridpoints  the  primitive  variables  u  and  not  the 
atixiliary  variables  w.  For  this  purpose  the  updates  in  w  and  (f>  corresponding 
to  the  relaxing  the  elliptic  factor  at  point  i  should  be  translated  into  the 
updates  in  u,p  at  points  {i  —  1),  (i),  {i  +  1)  according  to  M^. 


5.1.3  DGS  relaxation  and  the  Riemann  solver 

It  is  well  known  that  the  one-dimensional  Eider  system  can  be  diagonalized, 
i.e.  rewritten  as  a  set  of  (locally)  decoupled  advection  equations  for  the  char¬ 
acteristic  variables  (s,  a'^,  a  )^,  where 

Q!+  =  pcUx+Px,  and  =  pcUx  —  Px-  (49) 

Some  algebra  reveals  that  relaxing  the  Full-Potential  (elliptic)  factor  corre¬ 
sponds  to: 

-  at  point  (i  —  1)  -  update  Q!+,  keep  ot~  the  same; 

-  at  point  (i  -f  1)  -  update  cx~,  keep  O'"'"  the  same; 

-  at  point  i  -  update  both  O'"*'  and  a  . 

We  can  conclude  that  there  are  some  links  between  the  characteristic  variables 
formulation  (approximate  Riemann  solver)  and  the  design  of  the  DGS. 
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5.2  Two-dimensional  case 


We  shall  look  now  for  a  factorizable  upwind  scheme  in  two  dimensions. 


5.2.1  Dimension-by-dimension  approach 

Considering  here  the  isentropic  case  (no  loss  of  generality  for  the  purpose  of 
the  discussion  presented  here),  we  can  write  such  a  scheme  in  the  following 
symbolic  form 

=  0  (50) 


The  modified  equations  (or  the  First  Difierential  Approximation  -  FDA)  cor¬ 
responding  to  the  scheme  dimension-by-dimension 


FDA{L^)  = 

( p[~2  ®  ~  2  ^ 


0 


/3[-|(|«|4x  +  Cdyy)  -{■Q]dy-  \\vdyy 


(51) 


2c^yy^ 


-l-A  +  QJ 


It  is  easy  to  see  that  the  matrix  (51)  cannot  be  factorized. 


5.2.2  Genuinely  multidimensional  approach 

The  approach  towards  the  constniction  of  discrete  schemes  for  the  Euler  equa¬ 
tions  ([1,2])  was  called  “genuinely  multidimensional”  since  it  leads  to  schemes 
that  retain  (to  a  certain  extent)  the  rotational  invariance  property  of  the  Eu¬ 
ler  equations.  Namely,  it  can  be  argued  that  the  artificial  dissipation  terms 
present  in  these  schemes  approximate  a  rotationaUy  invariant  differential  oper¬ 
ator.  In  its  turn  this  may  mean  that  the  waves  oblique  to  the  grid  are  “properly 
upwinded”  or,  in  other  words,  the  same  approximate  Riemann  solver  can  be 
“recovered”  in  an  arbitrary  direction. 

It  is  not  clear  whether  or  not  this  property,  though  intuitively  appealing,  is 
of  any  direct  practical  importance  for  the  steady-state  computations.  There¬ 
fore,  we  do  not  discuss  it  in  detail.  However,  we  have  observed  previously  that 
there  are  some  links  between  the  approximate  Riemann  solver  and  the  design 
of  DGS  in  one  dimensional  case.  Therefore,  the  following  conjecture  seems 
reasonable. 
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Conjecture  A  genuinely  multidimensional  scheme  is  factorizable. 

It  was  pointed  out  in  [2]  that  some  of  the  multidimensional  second  order  correc¬ 
tions  can  be  added  without  limiters  resulting  in  a  linear  ‘positive’  scheme  with 
essential  multidimensional  character.  Namely,  for  the  u-momentum  equation 
in  subsonic  case,  those  are  the  cross-derivative  correction  terms  that  compen¬ 
sate  for  the  loss  of  accuracy  due  to  the  artificial  dissipation  in  x  direction.  For 
the  u-momentum  equations  those  will  be  the  terms  that  compensate  for  the 
loss  of  accuracy  due  to  the  artificial  dissipation  in  y-direction. 

Writing  such  a  scheme  in  the  symbolic  form 

4d]u"  =  0,  (52) 


and  considering  the  corresponding  FDA 


FDA{L%j,{)  = 

^  +  Q]  p[~\ip  ~  dx  — 

p\~\ip~  l'^l)^®y]  p[~\{\'^\^XX  +  cdyy)  +  Q]  dy—  2cQy 

\  p{dx  ~  2^9xx)  P{^y  ~  2c^yy) 


(53) 


we  can  easily  verified  that  FDA  matrix  is  factorizable.  The  added  multidimen¬ 
sional  correction  played  a  crucial  role  in  achieving  this  property. 


However,  the  FDA  does  not  uniquely  define  a  discrete  scheme.  Moreover,  not 

all  the  discrete  schemes  corresponding  to  a  certain  FDA  are  factorizable.  A 
factorizable  scheme  corresponding  the  above  mentioned  FDA  was  constructed 
on  a  quad-type  grid.  The  details  concerning  the  scheme  as  well  as  the  Dis¬ 
tributive  Gauss-Seidel  relaxation  will  be  given  elsewhere. 


6  NUMERICAL  EXPERIMENTS 

6. 1  Multidimensional  upwinding 


The  purpose  of  the  numerical  experiments  reported  in  this  section  is  to  demon¬ 
strate  the  robustness  of  the  genuinely  multidimensional  upwind  scheme  and 
the  quality  of  the  nmnerical  sohitions  obtained  by  its  means.  The  multigrid 
algorithm  employs  the  lexicographic  Gauss-Seidel  relaxation. 
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Fig.  3.  Supersonic  flow  in  a  channel  over  a  circular  bump,  grid  161  x  33:  a)  solution 
obtained  by  2FMF  -  W{2, 1)  algorithm;  b)  the  same,  except  that  3  more  W{2, 1) 
cycles  were  performed  on  the  finest  grid. 


It  was  demonstrated  in  [4], [5]  that  the  2FMG  -W{2,1)  algorithm  employing 
the  genuinely  two-dimensional  advection  .scheme  (and  Gauss-Seidel  relaxation 
with  direction-free  ordering  -  like  Red-Black)  is  capable  of  producing  second 
order  accurate  (both  in  smooth  regions  and  in  terms  of  discontinuity  location) 
sohition  to  stich  a  problem.  We  expect  this  to  be  true  for  the  Euler  equations 
as  well  (though  more  studies  .should  be  performed).  Therefore,  we  present 
solutions  obtained  u.sing  this  algorithm  for  a  few  testcases. 


6.1.1  Supersonic  flow  in  a  channel  with  a  bump. 

The  test  case  considered  here  is  a  supersonic  (Mach=2.9)  flow  in  channel  with 
a  circular  bump.  The  bump  is  located  at  the  lower  wall  of  the  channel  at 
1  <  a:  <  2  and  its  surface  is  a  circular  arch  of  r/S  and  radius  1.  Note  that 
the  actual  shape  of  the  domain  is  a  rectangle.  The  influence  of  the  bump  on 
the  flow  is  imposed  through  the  boundary  conditions:  the  velocity  component 
normal  to  the  surface  of  the  bump  at  a  certain  location  is  being  reflected. 

The  experiment  uses  the  finest  grid  of  the  size  161  x  33  points.  The  .solution 
obtained  by  2FMG  -  W(2, 1)  algorithm  is  pre.sented  on  Fig.3(a).  Fig.3(b) 
presents  the  numerical  solution  obtained  using  the  same  algorithm  but  per¬ 
forming  3  more  cycles  (total  5)  on  the  finest  level. 
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Fig.  4.  Transonic  flow  over  a  wall  with  a  circular  bump:  free  stream  Mach=.d,  grid 
200  X  200  pts. 

6.1.2  Transonic  flow  over  a  circular  bump 

The  testcase  considered  here  is  concerned  with  a  transonic  flow  over  a  flat 
wall  with  a  bump.  The  sxirfaee  of  the  bump  is  again  a  circular  arch  of  7r/3  and 
radius  1  and  its  location  is  between  3.5  <  a;  <  4.5.  Again,  in  order  to  keep  the 
experiments  simple  at  this  stage  of  work,  the  bump  is  treated  the  same  way 
as  in  the  previous  experiments.  The  free  stream  Mach=.9  in  this  case.  The 
shock  of  the  “fish-tail”  shape  is  generated  in  this  case  (Pig.4). 


6,1.3  Subcritical  flow  past  an  airfoil. 

5  Here  we  present  an  experiment  concerned  with  the  subcritical  fiow  past  an 
airfoil.  The  testcase  considered  is  Mach=  .63  flow  past  NACA0012  airfoil  at 
the  angle  of  attack  of  2°.  The  grid  conatins  abotit  9800  nodes.  Pressure  and 
density  contours  are  presented  by  Figs.5  (a)  and  (b)  respectively. 


6.2  Pressure  equation  based  schemes 


Here  we  illustrate  the  efficiency  of  the  multigrid  algorithm  that  distinguishes 
between  the  different  co-factors  of  the  system.  The  testcase  is  a  fiow  in  a 
channel  with  a  bump.  The  mutligrid  cycle  employed  is  P(2, 1),  the  relaxation 
is  lexicographic  Gauss- Seidel  in  the  flow  direction.  The  solution  contour  plots 
for  the  incompressible  case  are  presented  on  Fig.6.  The  sample  convergence 
rates  for  the  incompressible  and  compressible  subsonic  cases  can  be  found  in 
Fig.7.  It  can  be  obserevd  that  the  residuals  reduction  per  multigrid  cycle  is 
almost  an  order  of  magnitude.  Also,  the  extensive  set  of  computational  results 
presented  in  [21]  indicates  that  this  behavior  is  essentially  grid-independent. 
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a) 


b) 


Fig.  5.  Subcritical  flow  (Mach  =  .63)  over  NACA0012  aurfoil,  angle  of  attack  -  2°, 
grid  of  «  9800  nodes:  a)  pressure  contours;  b)  density  contours. 

7  CONCLUSIONS 


An  approach  towards  constructing  a  gemiinely  multidimensional  upwind  scheme 
was  intoduced  in  [1-3].  The  fundamental  advantage  of  this  approach  for  prac¬ 
tical  purposes  is  that  the  multidimensional  high-resolution  mechanism  (tmlike 
the  standard  one)  does  not  damage  the  stability  properties  of  the  scheme.  The 
conclusion  made  in  this  paper  is  that  the  genuinely  multidimensional  approach 
leads  to  a  scheme  that  is  also  factorizable.  The  practical  importance  of  this 
property  is  that  an  optimally  eflBcient  multigrid  solver  can  be  obtained  through 
the  constniction  of  an  appropriate  Distributive  Gauss-Seidel  relaxation,  that 
distinguishes  between  the  different  co-factors  of  the  equations.  Also,  since  the 
factorizability  property  can  be  easily  verified,  we  suggest  it  is  used  as  the 
definition  of  genuine  multidimensionality  of  a  scheme. 
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a) 


b) 


Fig.  6.  Incompressible  flow  in  a  channel  with  a  bump,  grid  of  97  x  33  vertices:  a) 
ii-velocity  contours;  b)  pressure  contours. 


Fig.  7.  Sample  convergence  histories  for  incompressible  and  compressible  (Mach=0.2 
and  Mach=0.5)  cases. 

We  also  presented  briefly  the  pressure  equation  based  schemes  and  demon¬ 
strated  on  their  example  what  efliciency  of  the  multigrid  solver  can  be  ex¬ 
pected  if  the  distinction  between  the  different  co-factors  of  the  equations  is 
made  by  the  algorithm. 

Due  to  its  generality,  the  genuinely  multidimensional  approach  for  discretiza¬ 
tion  of  the  Euler  equations  may  play  a  crucial  role  in  constructing  a  general 
optimally  efficient  multigrid  flow  solver  suitable  for  engineering  computations. 
This  is  because 

-  it  does  not  rely  on  casting  the  equations  into  any  special  form; 

-  extends  to  the  compressible  Navier-Stokes  equations. 
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